function out=f_segregation_lowF(F_centers,Combined_dual,ifplot,L_or_R,soft_hard,rad,f_cut)
PlotF = [8 20 40 80];
fs = 32.1e3; % this is the one for the 3D3A database
% rad = 0.3;% official 0.3
if ifplot==1
    f1=figure('position',[100 100 1500 500]);
    f2=figure('position',[100 100 400 250]);
    f3=figure('position',[500 400 400 250]);
end
x1=Combined_dual(:,1);
x2=Combined_dual(:,2);
nfft = round(40e-3*fs)
win = hamming(nfft,'periodic');
Combined_mono = sum(Combined_dual,2);
% y_ifft = spectrogram(Combined_mono);
% % % %      short time fourier transform:
% [SL,SR,freq_pos, flag01]=mySTFT(x1,x2,fs); % replace it with the MATLAB stft
[SL,Fmask,Tsig] = stft(x1,fs,'Window',win,'FFTLength',nfft);
if ifplot==1
    figure(f2)
    % subplot(131)
    stft(x1(1:fs*2),fs,'Window',win,'FFTLength',nfft);
    ylim([0 5]);title('A. Spectrogram of Target + Noise')
    colormap('jet')
    exportgraphics(gcf, 'AA_Figure3A_spectrog.png', 'Resolution', 300,'ContentType', 'vector')

end

% % % % There are 1284 (nfft=1284) frequencies in the stft, with 642 being
% nonredundant.
if ifplot==1
figure(f1)
end
SR = stft(x2,fs,'Window',win,'FFTLength',nfft);
Combined_stft = stft(Combined_mono,fs,'Window',win,'FFTLength',nfft);
size(Combined_stft)
% sdf0
freqs = (1:(nfft/2))./(nfft/2).*fs/2;
freqs(PlotF)
SL = SL(end/2+1:end,:);
SR = SR(end/2+1:end,:);
B_mask = ones(size(SL));
% 'size of bmask'
% size(B_mask)
% ser0
%      normalization:
% [SL1,SR1]=Normalize(SL,SR);
fcount = 0;
[SL1,SR1]=Normalize_Both(SL,SR,L_or_R);
for ifreql = 1:length(freqs)/2
    freq = freqs(ifreql);
    %     X_P = [real(SR1(ifreql,:))' imag(SR1(ifreql,:))'];
    if L_or_R ==2 % old
        X_P = [real(SR1(ifreql,:))' imag(SR1(ifreql,:))'];
    else
        X_P = [real(SL1(ifreql,:))' imag(SL1(ifreql,:))'];
    end
    for i_time=1:length(X_P)
        % B_mask(ifreql,i_time)=1;
        dist = sqrt((X_P(i_time,1)-F_centers(ifreql,1)).^2 + (X_P(i_time,2)-F_centers(ifreql,2)).^2);
        if dist < rad & freq>f_cut
            B_mask(ifreql,i_time)=soft_hard;
        end
    end
    if ifplot==1 & length(find(PlotF==ifreql))>0
        figure(f1)
        fcount=fcount+1
        subplot(1,4,fcount)

        Ind_M = find( B_mask(ifreql,:)==1);
        plot(X_P(Ind_M,1),X_P(Ind_M,2),'b.','Markersize',4);hold on;


        Ind_M = find( B_mask(ifreql,:)==0);
        % plot(X_P(Ind_M,1),X_P(Ind_M,2),'g.','Markersize',4);hold on;
        plot(X_P(Ind_M,1),X_P(Ind_M,2),'b.','Markersize',4);hold on;

        set(gca,'Xtick',-1:0.5:1)
        set(gca,'Ytick',-1:0.5:1)
        title([ num2str(freq), ' Hz'],'fontsize',14);
        xlabel('Real');ylabel('Imag')
        hold on;
        plot(-1:1,[0,0,0],'k', 'LineWidth', 0.5);
        plot([0,0,0],-1:1,'k','LineWidth',0.5)
        axis([-1.1 1.1 -1.1 1.1])
        axis square
        box off
        s=plot(F_centers(ifreql,1),F_centers(ifreql,2),'ro');set(s,'markersize',46,'linewidth',2)
        pause(0.01)
    end % plotting

end
full_mask = [flipud(B_mask); B_mask];
if ifplot==1
    figure(f3)
    contourf(Tsig,Fmask./1000,full_mask,'edgecolor','none');
    ylim([0 5]);
    xlabel('Time (s)');ylabel('Frequency (kHz)');colorbar
    caxis([0 1])
    colormap('gray')

    if soft_hard==0
        title('B. Hard Mask')
        exportgraphics(gcf, 'AA_Figure3B_hardmask.png', 'Resolution', 300,'ContentType', 'vector')
    else
        title('C. Soft Mask')
        exportgraphics(gcf, 'AA_Figure3C_softmask.png', 'Resolution', 300,'ContentType', 'vector')

    end
end
out.reconstructed = real(istft(Combined_stft.*full_mask,fs,'Window',win,'FFTLength',nfft));
out.freqs = freqs;
if ifplot==1
out.f1 = f1;
end
